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Abstract 

It is well-known that the dynamics of the Arnol'd circle map is phase-locked in regions of the 
parameter space called Arnol'd tongues. If the map is invertible, the only possible dynamics 
is either quasiperiodic motion, or phase-locked behavior with a unique attracting periodic 
orbit. Under the influence of quasiperiodic forcing the dynamics of the map changes dra- 
matically. Inside the Arnol'd tongues open regions of multistability exist, and the parameter 
dependency of the dynamics becomes rather complex. This paper discusses the bifurcation 
structure inside the Arnol'd tongue with zero rotation number and includes a study of nons- 
mooth bifurcations that happen for large nonlinearity in the region with strange nonchaotic 
attractors. 
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1 Introduction 



The Arnol'd circle map | Arnol'd, 1965 1 



x n+ i = x n + 9, + ^- sin(27rx„) mod 1 (1) 

is one of the paradigms for studying properties of nonlinear dynamical systems, both because 
it is a very simple map, and because of its great physical relevance (see, e.g. fBohr et al.,| 
1985[| ). Using the circle map one can model the structure of phase- lockings (devil's staircase) 



of a periodically forced nonlinear oscillator [ [Jensen et al., 1983 ; Jensen et al., 1984| and the 



current-voltage characteristics of a driven Josephson junction [ Bohr et al., 1984|]. The phase- 



locked regions of the Arnol'd circle map form the well-known Arnol'd tongues [ Arnol'd, 1983 



Hall, 1984]. If | if | < 1 there is a unique periodic attractor with a particular rotation number in 
each tongue. 

In this paper we study the structure of the phase-locked regions of the Arnol'd circle map 
driven by a rigid rotation with an irrational frequency. This system exhibits different kinds 
of dynamics, namely quasiperiodic motions with two and three incommensurate frequencies, 
chaotic attractors, and strange nonchaotic attractors (SNAs). SNAs have a strange geometrical 
structure, but unlike chaotic attractors they do not exhibit a sensitive dependence to changes 
in the initial conditions, i.e. their dynamics is not chaotic. They have been found in many 
quasiperiodically forced systems [|Grebogi et al., 1984 1, and also in the quasiperiodically forced 



circle map [Ping et aI.,198S| ; |Feudel et al., 1995(1 . 



Previous investigations show that regions of bistability occur in phase-locked regions of the 
quasiperiodically forced circle map [ Glendinning et al. , 2000(1 an d the phase-locked regions change 



in shape depending on the strength of the forcing. This change in shape is related to the emer- 



gence of SNAs [ Feudel et al., 1997 1- Based on these studies the aim of this paper is twofold. We 
study regions in parameter space where more than two attractors coexist (pockets of multista- 
bility). Secondly, we discuss the relation between these multistable regions and the appearance 
of strange nonchaotic attractors. In our discussion of the way these attractors are created and 
destroyed we are led to a description of nonstandard (nonsmooth) bifurcations of the invariant 
curves. 

If the unforced circle map is modified by introducing additional nonlinearities, coexisting 



attractors with the same rotation number can occur within the phase- locked regions [McGehee 



(fe Peckham, 19961 . In our system, multistability within the phase-locked regions is induced 



by the forcing rather than an additional nonlinear term. In quasiperiodically forced systems 
the coexisting attractors may be either invariant curves or SNAs depending on the strength of 
the forcing. Our investigation focuses on how these multistable regions appear and disappear 
under variation of the system's parameters: the nonlinearity K and the forcing amplitude e. 
For the tongue with zero rotation number the multistable regions open and close by smooth 
saddle-node or pitchfork bifurcations of invariant curves if K and e are small. For larger K 
these saddle-node and pitchfork bifurcations become nonsmooth: instead of merging uniformly 
(smooth bifurcation), the relevant stable and unstable invariant curves appear to collide only in 
a dense set of points. 

The paper is organized as follows. Section ^ recalls important properties of the Arnol'd circle 
map relevant for this study and their changes under the influence of quasiperiodic forcing. In 
particular we discuss the phase-locked region with zero rotation number. Within this phase- 
locked region we find pockets of multistability with a rather complex bifurcation structure which 
is analyzed in Sec. ||. Smooth and nonsmooth saddle-node and pitchfork bifurcations, leading to 
coexisting attractors, are studied in Sec. |I] to get a better understanding of the changes in the 
bifurcation structure depending on the strength of nonlinearity and forcing. Furthermore, we 
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investigate the transition between smooth and nonsmooth bifurcations and its implications to 
the dynamics of the system. In the full parameter space we find bifurcations of codimension two. 
In Sec. [| we discuss a special codimension-2 point that involves only nonsmooth bifurcations. 
It turns out that the unfolding of this point is very different from the smooth analog. Finally, 
in Sec. ^, we briefly discuss phase-locked regions with small, but finite, rotation number. We 
conclude this paper with a summary in Sec. [?] and an Appendix with details on the numerical 
computations. For readers with a black and white copy of this article we provide a supplementary 



website | Osinga et al., 200C | . 



2 The Quasiperiodically Forced Circle Map 

The quasiperiodically forced circle map is a map on the torus with lift 

K 

x n+1 = x n + O + — sin {2-Kx n ) + e sin (27n? n ) , (2) 

#n+i = ^n+w modi, (3) 

where i? n and x n modulo 1 give the coordinates on the torus. The parameter Q is the phase 
shift, K denotes the strength of nonlinearity (K > 0), e is the forcing amplitude, and the forcing 
frequency to is irrational. Throughout this paper we choose to work with u) = (y/5 — l)/2. 

2.1 The unforced system 

Let us recall the behavior of the unforced circle map (|T|). The dynamics of this map can be 
either periodic, quasiperiodic, or chaotic, depending on the parameters £1 and K. The critical 
line K = 1 divides the parameter space into two regions. If K < 1 the map is invertible and the 
motion can only be periodic (phase-locked) or quasiperiodic. For K > 1 the map is noninvertible 
and chaotic motion is possible. 

The rotation number is used to characterize the different kinds of motion. It is defined as 

p(n,K)= lim (4) 

N^oo iv 

where xn is the iVth iterate of (|T|), starting from xq. It can be shown that p(Q,K) does not 
depend on xq if K < 1. If the rotation number is rational, the attracting motion is periodic, 
otherwise it is quasiperiodic. For K < 1 the parameter space is split into regions with ratio- 
nal rotation number, the phase-locked regions or Arnol'd tongues, and regions with irrational 
rotation number corresponding to quasiperiodic motion. For example, the main tongue, the 
phase-locked region with zero rotation number, is bounded by the curves f2 = ±CIq(K), where 
Qq(K) = K/2tt. For any choice of Vt and K, with < f^o and K < 1, there are exactly two 
fixed points, one is attracting and the other is repelling. At the boundary \Q\ = f^o the two 
fixed points are annihilated in a saddle-node bifurcation. For other rotation numbers p ^ 
the iT-dependency of the boundary Q, p {K) is nonlinear, but it is always a curve of saddle- node 
bifurcations. 

2.2 The forced system 

A variety of behavior is possible in the coupled maps (|2])-(|3|). The rotation number (||) exists, 
but depends on e and uj in addition to ft and K, and the direct analogs of the periodic and 
quasiperiodic motion of the uncoupled Arnol'd map are invariant curves (the graph of a function 
$ i — > x(i))) and motion which is dense on the torus, respectively. In the former case the rotation 
number is rationally related to u (p = r\ + r2iv with rational, i = 1,2) and in the latter 
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Figure 1: The boundary — Qo(e,K) of the phase-locked region with zero rotation number; e € [0,5] runs 
from left to right, K £ [0,1] from back to front, and |fi| £ [0,0.16] from bottom to top. Pairs of invariant 
curves inside the phase- locked region (|fl| < flo) are annihilated at its boundary leading either to three- frequency 
quasiperiodic motion (yellow area) or SNAs (red area); see App. A.l. 



case there is no such rational relation. For both rationally and irrationally related rotation 
numbers strange nonchaotic attractors (SNAs) may also be possible. An SNA has a strange 
geometric structure, that is, it can be viewed as the graph of an everywhere discontinuous 
function $ i— ► but the dynamics on the attractor is not chaotic, because typical Lyapunov 

exponents in the x-direction are negative (there is always a zero Lyapunov exponent in the 
^-direction) . 

If the rotation number is rationally related to u then the motion is said to be phase-locked 
and the regions of parameter space in which the motion is phase-locked are analogous to the 
Arnol'd tongues of the unforced map. On the boundaries of the phase-locked regions we expect 



to see saddle- node bifurcations. There is an additional complication in the forced maps [|Feudel 
et al., 1995 1 in that the saddle-node bifurcations may be smooth (two invariant curves converge 
uniformly from inside the phase- locked region) or nonsmooth. In the smooth saddle-node bi- 
furcation the nontrivial Lyapunov exponent in the x— direction goes to zero at the bifurcation 
point. In the nonsmooth saddle-node bifurcation the two invariant curves appear to collide only 
on a dense set of points. Moreover, the typical nontrivial Lyapunov exponent remains negative. 
These nonsmooth saddle-node bifurcations seem to be associated with the appearance of SNAs 
outside the phase-locked region ||Ding et al.,T989| ; |Feudel et al., 1995| ; |Glendinning, 1998|] , 

It can be shown that one mechanism of the appearance of SNAs is related to changes in the 
shape of the phase- locked regions | Feudel et al., 1997 ; |Glendinning et al., 2000 1. For the unforced 
map the width of a phase-locked region increases monotonically with increasing nonlinearity K; 
this is no longer the case for positive forcing amplitude e. Moreover, as Fig. |l| shows, for 
fixed K the width of the phase-locked region oscillates as e increases. In particular, there are 
certain values of e for which the width of the phase-locked region becomes extremely small. 
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Figure 2: The bifurcation structure for K = 0.8 of the phase-locked region with zero rotation number (a) with 
detail of the first overlap (b). 



Unfortunately, using only numerical methods, we cannot decide whether the region actually 



closes or not at those e-values. For more details on numerical computations we refer to App. |A.1| . 

For small fixed nonlinearity K the boundary £Iq{e, K) of the phase-locked region can be 
approximated by (the modulus of) a Bessel function of order zero using first order perturbation 



theory [ Glendinning et al., 2000] ] . Numerical simulations also revealed regions of bistability in 



the vicinity of the zeroes of the Bessel function, where the width of the phase-locked region is 
very small. The bistability regions are bounded by saddle-node bifurcations of invariant curves, 



which has been confirmed by second order perturbation theory Glendinning fc Wiersig, 1999 ]. 



The study in | Glendinning et al., 2000| ; Glendinning Sz Wiersig, 1999| only applies for K close 



to 0. We wish to study what happens to the phase-locked region with zero rotation number 
for larger K. However, we restrict our considerations to the invertible case K < 1, so that 
chaos is ruled out. We find that the regions of bistability contain other regions where even more 
attractors coexist. In the following we describe how these regions appear and disappear as a 
parameter varies. We also study smooth and nonsmooth bifurcations and make some remarks 
on the appearance of SNAs. 

The majority of the rest of this paper describes the results of numerical simulations of 
Eqs. As such, the reader should bear in mind that our conclusions are based on numerical 

observations and may turn out to be misleading in places. We have made every effort to avoid 
such problems (see the Appendix) and believe that the phenomena reported are sufficiently 
interesting and mathematically intractable to merit this numerical investigation, even if we 
remain uncertain of some of the outcomes. The reader is encouraged to maintain a healthy 
scepticism throughout. 

3 The Internal Structure of the Main Tongue for K = 0.8 

In the simplest case the boundary \Q\ = £Iq(e,K) represents the disappearance of two invariant 
curves, a stable and an unstable one. However, for |f2| < Qq(e,K) more than two invariant 
curves may exist that disappear before this boundary is crossed. Such pockets of multistability 
are found near local minima of ^Iq(e,K), cf. the region of bistability predicted by perturbation 
analysis for small K [ Glendinning et al., 200C ; |Glendinning fc Wiersig, 1999 ] . For example, 



Fig. |2|(a) shows a cross-section of Fig. |T| at K = 0.8, with both positive and negative sides of the 
boundary of the phase-locked region. The outer boundary is the function \D\ = 0,q(e, 0.8). Extra 
curves are drawn marking the boundaries of pockets of multistability, which is best seen in the 
enlargements Figs. §(b), |] and ||. These pockets of multistability can be considered as overlaps 
of different "bubbles" with the same rotation number as in Fig. @(a). For better visualization 
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(a) (b) 

Figure 3: Sketch of the bifurcations along the line = and K = 0.8 in the first (a) and second (b) overlap. 
Shown are invariant circles (represented by one point) versus e. Closed curves represent stable and dashed curves 
represent unstable circles. The colors correspond with the colors of the bifurcation curves in Figs. || and | 




Figure 4: The bifurcation structure for K = 0.8 in the second overlap (a) seems to be the same as in Fig. Q(b). 
However, an enlargement (b) shows that the structure is much more complicated. 



we have chosen different colors for bifurcations of different pairs of invariant curves. 

3.1 Bifurcations for K = 0.8 in the first region of overlap 

The first overlap is enlarged in Fig. |2](b). The orange and purple curves enclose a rhombus 
shaped region where two attracting and two repelling invariant circles exist. This region of 
bistability is bounded by curves of saddle-node bifurcations that end in pitchfork bifurcations 
on the line f2 = 0. Note that if O = the map has a symmetry {x \— * —x, i? i— ► $ + 1/2) which 
implies that the rotation number in the x-direction is always zero and pitchfork bifurcations 
should be expected. The bifurcation sequence for Q = is sketched in Fig. ||(a) where each 
circle is represented as a point and e increases along the horizontal axis. The bottom and top 
lines are identical, representing the modulo 1 computations. The purple and orange dots are 
the pitchfork bifurcations that mark the crossing of the purple and orange curves in Fig. ||(b) 
along £2 = 0, respectively. 

3.2 Bifurcations for K = 0.8 in the second region of overlap 

Figure |I](a) shows a detail of the second overlap. This picture is very similar to Fig. ||(b), but 
the bifurcation diagram along 17 = in Fig. |3|(b) reveals a more complex structure. Let us first 
discuss Fig. ||(b). The first bifurcation (light blue dot) is the same pitchfork bifurcation as the 
purple dot in Fig. ||(a): the stable circle becomes unstable, creating two new stable circles. As 
e increases two other stable circles are born in a pair of saddle- node bifurcations (red dots). We 



MULTISTABILITY AND NONSMOOTH BIFURCATIONS 



7 




now have four different attractors. Note that these saddle-node bifurcations happen at the same 
values of e due to the symmetry along the line f2 = 0, as referred to earlier. 

At the purple dots the two attracting circles from the pitchfork bifurcation disappear in a 
pair of saddle-node bifurcations. Note that this pair of saddle-node bifurcations is connected to 
the purple pitchfork bifurcation of Fig. |3|(a) via the purple curve off Q = in Fig. ||(a). We are 
now left with two attractors and two repellors. These last two attractors do not disappear in a 
pitchfork bifurcation as in Fig. ||(a) (orange dot). Instead, they disappear in a pair of saddle- 
node bifurcations (dark blue dots) with two repellors that are born in a pitchfork bifurcation 
(green dot) for slightly smaller e. Note that this pitchfork bifurcation is subcritical, as opposed 
to the supercritical orange one in Fig. ||(a). 

The unfolding of these bifurcations inside the phase-locked region with £1 ^ is shown in 
detail in Figs. |] and [|. The curves are colored according to the colors of the bifurcations in 
Fig. ||(b). As expected, for Q ^ pairs of saddle-node bifurcations no longer happen at the same 
values of e. They form two different curves that cross each other exactly at f2 = 0. We already 
mentioned earlier that the purple curves connect all the way left in a pitchfork bifurcation on 
the line f2 = in the first overlap. The light blue curves start in the pitchfork bifurcation at 
O = and become the outer boundary |fi| = $7o( £ >0.8) once they cross the purple curves. The 
red curves form swallowtails with the purple curves on the right side and the dark blue curves 
on the left side; see Fig. 0(b). Finally, the dark blue curves form swallowtails with red and green 
curves; see Fig. ||. 



4 The Structure of Bifurcations in the (e, K)-plane 

The bifurcation structure depends on the strength of the nonlinearity K. We study the two- 
parameter dependence only on the cross-section f2 = 0, because the unfolding inside the phase- 
locked region with $7 ^ is similar to that discussed in Sec. ||. As expected, the regions of 
overlap change shape with K. In particular, only the regions with no more than two attractors 
persist for small K and moderate e. This is shown in Fig. ^with a cross-section at O = in the 



(e, iC)-plane of the second overlap; see App. A. 2 for details on how this picture was generated. 



In the following sections we describe the bifurcations in more detail. Section 4.1 discusses 
the sequence of bifurcations that happen as K decreases. For small K (less than approximately 
0.8 for the bifurcations we have looked at), saddle-node and pitchfork bifurcations happen via a 
uniform collision of invariant curves: at the moment of bifurcation, two (saddle-node) or three 
(pitchfork) curves merge at each value of •&. We call these bifurcations smooth bifurcations. 
For K close to 1 the attractors may become extremely wrinkled, which gives rise to nonsmooth 
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0.6 ■ 



0.4 ■ 

large phase sensitivity 
4 invariant curves 
n -2 6 invariant curves 

■ 8 invariant curves 

I — — — — •— — 
1.555 1.575 1.595 1.615 1.635 1.655 



Figure 6: A cross-section at fl = of the second overlap over a 200 x 200 grid. A point in the white region 
corresponds to two invariant curves, one is stable and t he ot her is unstable. The region in which the stable curve 
is extremely wrinkled (large phase sensitivity, see App. A. 2) is marked as well. 
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(a) (b) (c) 

Figure 7: Sketches of the bifurcations along the line Q — in the second overlapping region as K decreases from 
0.8 to 0. The bifurcation structure of Fig. g(b) for K = 0.8 transforms into one like Fig. ^(a) by "absorbing" 
pairs of saddle-node bifurcations in the second pitchfork bifurcation. 



bifurcations: at the moment of bifurcation the invariant curves now collide only in a dense set 
of -^-values. The nonsmooth pitchfork and nonsmooth saddle-node bifurcations are described in 
detail in Sees. 4.2 and 4.3, respectively. 



4.1 Smooth bifurcations for f2 = in the second region of overlap 

Figure |6| indicates that at most two attractors exist for small K and moderate e. This means 
that for fixed small K, the bifurcation portrait looks like Fig. |3|(a). Hence, as we decrease K 
from K = 0.8 to 0, the extra pairs of saddle-node bifurcations (see Fig. ||(b)) need to disappear 
somehow. It turns out that the last pitchfork bifurcation, the green dot in Fig. ||(b), "absorbs" 
these saddle-node bifurcations one by one. In doing so, the pitchfork bifurcation switches from 
subcritical to supercritical and vice versa (a standard codimension-2 bifurcation). A sketch of 
this process along the line O = is shown in Fig. @(a)-(c). 

In the (0,e)-plane the picture changes as follows. The first swap, Fig. 0(a), comes about 
as the pitchfork point and the intersection point of the dark blue curves on O = collapse; see 
Fig. ||(b). When K decreases, these points move closer together, causing the slope of the dark 
blue curves to become steeper and the ends of the swallowtail to move closer to = 0. Upon 
collision the green curves disappear and the dark blue curves end in a supercritical pitchfork 
bifurcation. 

In the second swap the dark blue curves disappear in a similar way via a collision of the 
pitchfork point and the intersection point of the two red curves, making the pitchfork subcritical 
again. Note that, in order for this to happen, the pitchfork point crosses the intersection point of 
the two purple curves; compare Figs. ||(b) and ||(a). Figures ^(a)-(b) show why this is a crossing 
and not a collision: the pair of purple saddle-node bifurcations happens "far out" in state space 
from the pitchfork bifurcation. Therefore, the crossing is only a crossing in this projection on 
the (0,e)-plane. 

The third swap, Fig. 0(c), is identical to the first, causing the disappearance of the red 
curves. In this bifurcation diagram at most two attractors coexist, which is the desired situation 
for K small. 



4.2 Nonsmooth pitchfork bifurcations 

For K close to 1 the situation is more complicated, because some of the invariant curves are very 
wrinkled and the pitchfork bifurcation becomes nonsmooth. This bifurcation has been found by 
Sturman [1999] in a similar map. Let us discuss what happens for $7 = 0, along the line K = 0.9 
as we approach the pitchfork bifurcation by decreasing e, starting in the yellow region in Fig. ||. 
Figure |8|(a) shows all invariant curves just before the bifurcation. The two stable invariant 
curves (black and blue) correspond to the two outer branches of the pitchfork, the unstable 
invariant curve (red) relates to the inner branch of the pitchfork separating the two outer ones. 
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(b) 



Figure 8: Nonsmooth pitchfork bifurcation with fi = and K = 0.9: (a) before the collision (e = 1.56765): 
attractors in blue and black, unstable invariant curves in green and red; (b) after the collision (e = 1.5675): SNA 
(red) and the unstable invariant curve not taking part in the bifurcation (green). 



The fourth invariant curve (green) is also unstable, but it is "far away" and does not take part in 
the bifurcation. As we decrease e towards the bifurcation point the three invariant curves (blue, 
black and red) approach each other, but due to their wrinkled structures they appear to collide 
only in a dense set of ^-values instead of merging uniformly as in a smooth bifurcation. This 
indicates that at the moment of bifurcation the attractor is an SNA. Numerical evidence suggests 
that this SNA persists and smoothes out to an invariant curve over a small e-interval; see for 
example the red attractor in Fig. B(b) for e below the bifurcation value. It is possible that we 



see the reverse of fractalization, a mechanism for the appearance of SNAs reported in [ Nishikawa 
& Kaneko, 1996]. However, if we use the method of rational approximations for testing whether 



the attractor is an SNA, we get conflicting results; see App. [A.3 for more details. We remark 
that we get these conflicting numerical results only for the nonsmooth pitchfork bifurcation. In 
any case, after a further decrease in e the attractor is clearly a smooth invariant curve. 

It is important to note that the nonsmooth pitchfork bifurcation is uniquely defined as the 
moment of collision of three invariant curves and the locus of bifurcation lies on a curve in the 
(e, iQ-plane. The process of fractalization is a gradual process where the moment of transition 
from an invariant curve to an SNA is not well-defined numerically. We wish to emphasize that 
it is, therefore, completely unclear whether the set of parameter pairs (e, K) with £1 = that 
exhibit SNAs after the nonsmooth pitchfork bifurcation has zero or finite size. 

Since the boundary of the region of bi- or multistability for O = is given by smooth 
and nonsmooth pitchfork bifurcations there is a codimension-2 point in the (e, K )-plane where 
the smooth and the nonsmooth pitchfork bifurcation curves meet. An approximation of this 
codimension-2 point is (e, K) = (1.564, 0.89). In any neighborhood of this point we always find 
all three kinds of dynamical behaviors: one stable invariant curve, two stable invariant curves, 
and one SNA; although the latter may only exist on the nonsmooth bifurcation curve itself. 
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4.3 Nonsmooth saddle-node bifurcations 

If K is small and e moderate then the saddle-node bifurcations observed numerically involve 
two invariant curves on the torus which converge and destroy each other. At larger values of K, 
simulations suggest that two invariant curves touch on an orbit at the bifurcation point, so points 
of intersection of these sets are dense on the curves. (Strictly speaking, the invariant sets are no 
longer continuous at the bifurcation point, but we will continue to refer to them as curves.) For 
quasiperiodically forced circle maps we can distinguish two types of these nonsmooth saddle- 
node bifurcations: one-sided and two-sided. In the one-sided nonsmooth saddle-node bifurcation, 
the collisions occur between pairs of invariant curves on the cylinder. An example is shown in 
Fig. ||(a) where the invariant curves and some of their translates by one in the x— direction are 
computed close to the bifurcation point. In the two-sided case each stable invariant curve on 
the cylinder touches both the unstable invariant curve immediately above it and the unstable 
invariant curve immediately below it. On the torus this implies that at the bifurcation point 
the attractor is everywhere discontinuous. An example of a two-sided nonsmooth saddle-node 
bifurcation is shown in Fig. |lO|(a ). These two-sided nonsmooth saddle-node bifurcations are 



described in |Feudel et al., 1995 1, where it is shown that after the bifurcation (with 17 = 
and K fixed) the map has an SNA with unbounded motion in the x— direction (Fig. |i~0|(b)) 
despite the fact that the rotation number remains zero due to the symmetry of Eqs. ©-(jH) 



when $7 = 0. This implies that the diffusion in the x— direction is extremely slow; see [|Feude] 



et al., 1995 ; Stark et al., 1999 1 for further details. In general, the two-sided nonsmooth saddle- 



node bifurcation is of codimension two, but it occurs as a codimension-1 phenomenon due to 



the symmetry if $7 = 0. Figure 11 shows the range of dynamics observed in the third region of 
overlap in the plane Vt = 0. Unbounded SNAs are observed in the blue regions and two-sided 
nonsmooth saddle-node bifurcations occur on the boundary between the white and blue regions. 

The unbounded SNA of Fig. [H](b) must contain orbits which are unbounded above and 
orbits which are unbounded below |Stark et al., 1999 1. This bidirectional diffusive motion of the 



unbounded SNA with 17 = becomes effectively unidirectional if |17| is very small, leading to a 
nonzero rotation number | Stark et al., 199"9| . This suggests that these unbounded SNAs lie on 



the boundary of the phase-locked region. Fig. [T^ shows this boundary in (e, K, 17)- space. It is 
clear that the height (i.e. the width in 17) of the boundary is very small, if not zero, in regions 
of the (e, K) plane with 17 = which have unbounded SNAs (compare the low plateau on the 
left of Fig. [l^ with the region of unbounded SNAs of Fig. [ll|). As the height of the boundary 
becomes non-negligible we observe that the saddle-node bifurcation on the boundary has become 
one-sided (see Fig. |9|), and we believe that it remains one-sided and nonsmooth throughout the 



red areas of Fig. 12 with non-negligible height. These red regions of the boundary correspond to 



saddle-node bifurcations with negative nontrivial Lyapunov exponents, and appear to preceed 



the creation of SNAs with nonzero rotation numbers outside the phase- locked region [Ding et al., 
1989| ; |Glendinning, 199g] . 



If the phase-locked region really has zero height on the plateau, then it is wrong to refer 
to saddle-node bifurcations on the interior of the plateau: these points would correspond to 
a transition from an SNA with negative rotation number to an SNA with positive rotation 
number through an SNA with zero rotation number as 17 increases through zero. In this full 
three-parameter unfolding the two-sided nonsmooth saddle-node bifurcations are of codimension 
two, occurring on curves bounding the plateau and separating parts of the boundary of the 
phase-locked region with unbounded SNAs from parts with one-sided nonsmooth saddle-node 
bifurcations. 
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0.25 



Figure 9: One-sided nonsmooth saddle-node bifurcation with Q — 0.001 and e = 2.58: (a) before the collision 
(K = 0.865): attractors in black, blue and green, repellors in red, purple and orange; (b) after the collision 
[K — 0.868): attractor with the same initial condition as the black attracting invariant curve in (a); repellor with 
the same initial condition as the purple repelling invariant curve. The attractor moves upwards in forward time 
and the repellor moves downwards in reverse time. The rotation number is nonzero. 
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Figure 10: Two-sided nonsmooth saddle-node bifurcation with = and e = 2.58 (cf. the transition from the 
white to the blue region in Fig. |ll]): (a) before the collision (K — 0.927): attractors in black, blue and green, 
repellors in red, purple and orange; (b) after the collision (K = 0.928): SNA with the same initial condition as 
the black attracting invariant curve in (a). The strange repellor is not shown. 
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Figure 11: Section Q, = of the third overlap region. In the region of large phase sensitivity the attractor is 
very wrinkled or even an SNA. 
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Figure 12: A small part of the boundary = Qo(e,K) of the phase-locked region with zero rotation number; 
e G [2.52, 2.62]_runs from right to left, K G [0, 1] also from right to left, and |fi| G [0, 0.004] from bottom to top; 
compare Fig. A grid of 80 x 160 points in the (e, AT)-plane is taken. Red regions with negative Lyapunov 
exponent (using A < —0.005 as criterion) correspond to transitions to SNAs (the yellow spots are due to the finite 
threshold) . 
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(a) (b) 



Figure 13: (a) Magnification of Fig. [TTJ. The white region corresponds to pairs of invariant curves, regardless of 
their phase sensitivity, (b) Magnification of (a). 



5 Nonsmooth bifurcation points of codimension two 

In the full parameter space, and also already in the (e, iT)-plane with f2 = 0, we expect 
codimension-2 bifurcations. We have already seen some of these; for example, the transition 
from yellow to red in Figs. |l] and 12 is the codimension-2 bifurcation curve marking the transi- 



tion from a smooth saddle-node bifurcation to a nonsmooth one, respectively [Kuznetsov et al., 



200C 1 . Another example is the codimension-2 point, mentioned in Sec. 4.2 , where the smooth 
and nonsmooth pitchfork bifurcation curves meet. Note that there is no curve of codimension-2 
points in this case, since the pitchfork bifurcation is restricted to the plane O = 0. As discussed 
in the previous section, there is also the two-sided nonsmooth saddle-node bifurcation curve. 

In this section we want to draw attention to an interesting codimension-2 point in the plane 
$7 = 0. This point can be seen in Fig. |ll] and in the enlargements Figs. |l3|(a)-(b) as the 
point where the region with unbounded SNAs (blue) and the bistable region (yellow) touch. It 
can be characterized as the moment where two nonsmooth (supercritical) pitchfork bifurcations 
happen simultaneously, i.e. a region of overlap is pulled apart; compare also the sketch in Fig. 7 
of flGlendinning et al., 2000] . 

5.1 Smooth analog of the nonsmooth codimension-2 bifurcation point 

Let us first discuss the smooth analog of this nonsmooth codimension-2 bifurcation point. Sup- 
pose for f2 = and some K < 1 fixed the bifurcation diagram involves only two supercritical 
pitchfork bifurcations that occur in the order as shown in Fig. [11(a). Now assume that as we 
increase K, this order is switched before we reach K = 1, without changing the type of pitch- 
fork bifurcation from supercritical to subcritical. As shown in Fig. 0(b), this means that we 
necessarily need extra curves of saddle-node bifurcations. 

The smooth analog of the nonsmooth codimension-2 bifurcation point is the point where 
the two pitchfork bifurcations happen at the same parameter values. In the (e, iT)-plane the 



complete bifurcation diagram should look like Fig. 15. The two supercritical pitchfork bifurcation 
curves are colored light-blue and magenta. The green and light-green curves are saddle-node 
bifurcation curves; compare also the colors in Figs. |l4|(a)-(b). Above the green and light-green 
curves, but below the light-blue and magenta curves there are four attractors. If we cross either 
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(a) (b) 

Figure 14: Two different bifurcation diagrams along the line = that both involve two (smooth) supercritical 
pitchfork bifurcations. Shown are invariant curves (represented by one point) versus e with K fixed. If the 
pitchfork bifurcation creating two attractors does not happen first, we necessarily need extra curves of saddle- 
node bifurcations (b). 



the light-blue or the magenta curve from this region there are three attractors. Above the light- 
blue and magenta curves there are two attractors. We distinguish five qualitatively different 
K- intervals numbered I-V in Fig. 15. Interval I corresponds to Fig. 14|(a) and V to Fig. 0(b). 



The qualitative behavior in the intervals II— IV is given in Figs. [l£(a)-(c), respectively. The 
codimension-2 point that we are discussing here is the intersection point of the light-blue and 
magenta pitchfork bifurcation curves. The intersections of the green and light-blue curves, and 
the light-green and magenta curves are only intersections in this projection onto the (e, X)-plane 
as can be seen in Fig. |l^, where the transition from Fig. |l~6|(a) to (b) marks the crossing of green 
and light-blue, and the transition from Fig. |l6|(b) to (c) represents the crossing of light-green 
and magenta. 

5.2 A nonsmooth bifurcation point of codimension two 

The nonsmooth version of Fig. ^ looks surprisingly simple in contrast; see Fig. 17. The coloring 



of the bifurcation curves is as in Fig. 15 with the restriction that all curves represent nonsmooth 
bifurcations. The upper pair of bifurcation curves corresponds to two-sided nonsmooth saddle- 
node bifurcations and the lower pair to nonsmooth pitchfork bifurcations. 

We now have four regions with distinctive dynamics. To the left and to the right of the 
two pitchfork bifurcation curves there is one attractor. In between the two curves of pitchfork 
bifurcations there are two attractors. This is illustrated for a particular choice of the parameters 
in Fig. |l~8|(b). The two attractors (black and blue) are two invariant curves separated by unstable 
invariant curves (green and red). In contrast to the previously discussed case of the nonsmooth 
pitchfork bifurcation, both unstable invariant curves are now close to the attracting invariant 
curves, because we are near the codimension-2 bifurcation point. The two different pitchfork 



bifurcations are illustrated in Figs. |l8|(a) and (c). If we decrease the forcing amplitude e, the 
black, blue and green invariant curves disappear in a nonsmooth pitchfork bifurcation to form 
the green attractor coexisting with the red unstable invariant curve, as shown in Fig. |l8|(a). This 
corresponds to a crossing of the light-blue nonsmooth pitchfork bifurcation curve in Fig. [17]. We 
cross the magenta nonsmooth pitchfork bifurcation curve in this figure by increasing e. In this 
case the red, black and blue invariant curves disappear and form the red attractor coexisting with 
the green unstable invariant curve as in Fig. |l8|(c). The fourth region between the nonsmooth 
saddle- node bifurcations is characterized by the existence of unbounded SNAs. 

Finally, we remark that the SNA region seems to have a fractal-like structure in the neighbor- 
hood of the codimension-2 point; compare Figs. |l3|(a)-(b). This would imply that the boundary 



MULTISTABILITY AND NONSMOOTH BIFURCATIONS 



17 




(a) (b) (c) 

Figure 16: Qualitative bifurcation diagrams for A'- values in the K- intervals II (a), III (b), and IV (c) in Fig. [l^. 
The equivalent diagrams for the 7f -intervals I and V are shown in Fig. Fl4l 



of the phase-locked region is fractal! 



6 Tongues with Nonzero Rotation Numbers 

We studied the structure of the phase-locked region with zero rotation number in great detail. 
In this section we briefly discuss the geometry of the tongues with nonzero rotation numbers 
that are close to the main tongue. More precisely, we determine the boundaries of the tongues 
with rotation numbers p = 1/Ff., with the kth Fibonacci number = Fk_\ + i^_2 and 



F\ = F% = 1, using the same numerical procedure as before; see App. A^l. Several tongues for 
K = 0.99 are shown in Fig. [K]. The fluctuations of the widths of the tongues with p > are 
due to numerical errors, which are of the same magnitude as the widths themselves. The thick 
borderline £lo(e,K) in this parameter regime corresponds to transitions to SNAs or bifurcations 
of SNAs of the type discussed in Sec. [D| 

How these tongues approach the boundary Qq(e,K) for fixed (e, K) as k goes to infinity 
depends in general on e and K; compare e = 2.55 and e = 2.61 in Fig. This behavior can be 
quantified with a scaling law for the distance Q p — Qq between the main tongue and the tongues 
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£ 



Figure 17: Unfolding of the nonsmooth co dimension- 2 bifurcation point. The green and light-green curves cor- 
respond to nonsmooth saddle-node bifurcations, while the light-blue and the magenta curves represent nonsmooth 
pitchfork bifurcations. 



with nonzero rotation numbers p. We used several parameter pairs (e, K) in Fig. 2C and our 
numerical calculations strongly suggest the scaling law p ~ £l p for e and K such that S7o = 0, 
and p ~ y/Qp — Qq otherwise. This agrees with the linear scaling in the high-e limit conjectured 
by Ding et al. [1989], the square-root scaling in the Arnol'd circle map for e = and K > (see 
e.g. [MacKay & Trcsser, 1984]) and the trivial linear scaling of the pure rotation for £ = K = 0. 
The square root scaling is associated with the one-sided saddle-node bifurcation, and the linear 
scaling is related to the change in the rotation number of the SNAs as in Sec. 4.3. The change 
between these scalings is associated with a two-sided nonsmooth saddle-node bifurcation point 
of codimension two at which two one-sided nonsmooth saddle-node bifurcation curves meet 
(presumably the nonsmooth analog of a cusp bifurcation point). A detailed study of this curve 
of codimension-2 points could be of interest. 



7 Summary 

We have studied the structure of the phase-locked regions in the quasiperiodically forced circle 
map. In particular, we have found regions of multistability where several attractors coexist. 
These regions of multistability appear due to the emergence of additional pairs of invariant 
curves as a result of saddle-node or pitchfork bifurcations under the variation of the forcing 
amplitude. As a result, these regions look like overlaps of phase-locked regions with the same 
rotation number. 

Opening and closing of these pockets of multistability are due to saddle-node and pitch- 
fork bifurcations of invariant curves. These bifurcations can be either smooth or nonsmooth 
depending on the strength of nonlinearity and the forcing amplitude. This is organized by the 
type of interaction between the stable and unstable invariant curves. In the smooth case these 
curves approach each other uniformly and then touch uniformly in each value of i9, so that the 
bifurcation looks like a simple merging of the invariant curves, analogous to the unforced case. 
Nonsmooth saddle-node bifurcations appear due to a wrinkled structure of the participating 
invariant curves, which collide at the bifurcation only in a dense set of ^-values. The result of 
this bifurcation is the emergence of a strange nonchaotic attractor. Similar to the nonsmooth 
saddle- node bifurcation we find a nonsmooth pitchfork bifurcation, but the details of this bi- 
furcation are still unclear. Both for the saddle-node and the pitchfork bifurcation there are 
codimension-2 points in parameter space marking the transition from a smooth to a nonsmooth 
bifurcation. The exact determination of the codimension-2 pitchfork bifurcation point and the 
self-similarity properties in its neighborhood should be possible by using renormalization group 
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Figure 18: Nonsmooth pitchfork bifurcations in the neighborhood of the codimension-2 point K — 0.83811, 
tt = and e = 2.55483: (a) to the left of the bistable region (K = 0.838, e = 2.5542); (b) within the bistable 
region close to the higher codimension point (K = 0.838, e = 2.55483); (c) to the right of the bistable region 
(K = 0.838, £ = 2.5551). 



0.008 



0.006 




8 



Fi gure 19; Section K = 0.99 of the r egio ns with rotation numbers p = 1/Fk, k being 11, . . . , 16, close to the 
(grey) region with p = 0; compare Figs. |l lHjj. 



20 



OSINGA, WlERSIG, GLENDINNING & FEUDEL 



Figure 20: The distance Q p - fi ( 



-5 



Q. 




-7 
-8 

-13 -11 -9 

In (£2-Q ) 

of the tongue with rotation number p 



l/Fk to the main tongue, where k 



ranges from 11 to 18; compare also Fig. [19[ The data can be approximated by straight lines in a In — In plot with 
slopes 1.032 ± 0.013 for e = 2.55 (squares) and 0.488 ± 0.025 for e = 2.61 (stars). The nonlinearity K is in both 
cases equal to 0.99. 



techniques. 

We have also investigated a nonsmooth codimension-2 bifurcation involving the merging of 
two nonsmooth pitchfork bifurcations and two nonsmooth saddle-node bifurcations which leads 
to regions of unbounded SNAs. This has no straightforward analog in smooth bifurcations. 

The positions of phase-locked regions in the neighborhood of the region with zero rotation 
number were also described. Fixing the nonlinearity K and the forcing e, we found that the 
rotation number p of these regions scales linearly with Q p whenever the width of the main tongue 
appears to be zero, and as — Qq otherwise. 
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Appendix: Details on the Numerical Computations 

We now explain our numerical computations in more detail and give values for the accuracy 
parameters. The Appendix is organized such that each section is related to one section in the 



main text. Appendix A.l explains how to compute the tongue boundary for the main tongue 
with zero rotation number and relates to Sec. |2.2| . Appendix [A.2| discusses the generation of 
Fig. H in Sec. ||[ Here, we also discuss how to identify SNAs using the phase sensitivity exponent, 
which is related to the derivative of x n with respect to the external phase #. The numerical 
issues that arise when determining the nonsmooth pitchfork bifurcation and the accompanying 
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fractalization process, as reported in Sec. 4.2 , are described in more detail in App. [A3| , In 
contrast to the method for identifying SNAs in the second region of overlap as described in 



App. A. 2, we used another more efficient method for determining the regions where SNAs exist 



in the third region of overlap. This method is explained in App. A. 4. 



A.l Numerical computation of the boundary of the phase- locked region 

The boundary f2 = £lo(e,K) of the phase-locked region with p = is half its width due to the 
symmetry. It is approximated by estimating the boundary point 0,q on a grid of 320 x 40 points 
in the (e, K)-p\ane. Following [[Stark et al., 1999 1 we determine the rotation number @) within 
an accuracy of ±X/N by averaging over a sample of 25 orbits of length N = F28 = 317811 (after 
1000 preiterations to eliminate the effect of transients), where are the Fibonacci numbers 
F\ = F2 = 1 and Fk = -Ffc-i + Fk-2- The Fibonacci numbers Fk are used since the value 
•&F k after Fk iterations is close to the initial value i?o due to the fact that ratios of Fibonacci 
numbers are good rational approximants of our irrational driving frequency lo. The initial 
fwl = [0,0.2] is repeatedly bisected, preserving the relation p(fig ) < X/N < p(£Iq) 
until Qq — Qq < A$7 = 10~ 5 . Finally, we choose to be the 
mean value of Qq and Qq, or zero if the mean value is smaller than our numerical accuracy Afi. 



interval [f2 , * o Q 
to ensure that Qq E [Qq , 



Note that, as remarked in the Sec. 2.2, the question of whether S7o really vanishes or is just very 
small cannot be answered by using only numerical methods. 

To distinguish between smooth and nonsmooth saddle-node bifurcations of invariant curves 
on the boundary of the phase-locked region we compute the nontrivial Lyapunov exponent 



N-l 

X(Sl,e,K) = lim — > In 



n=0 



dx n 



+ 1 



8Xr 



lim — ^ In 1 1 + K 



7V->oo 



cos 2ttx t 



(5) 



n=0 



Vanishing A (yellow regions in Fig. [l]) indicates smooth saddle-node bifurcation while negative A 
indicates nonsmooth saddle-node bifurcation (red regions); see Feudel et al. [1997] and Sec. 2.2 . 



A. 2 The cross-section = of the second region of overlap 

Figure || shows the bifurcation structure at a cross-section (7 = 0. The picture was generated as 
follows. For each grid point N = F32 = 2 178 309 iterations of Eqs. (|2|)-(|3|) were computed using 
25 different initial conditions (xq,&o = 0). We take advantage of the fact that inside the tongue 
with zero rotation number, each attractor is represented as a single- valued function x = X($), 
1? G [0,1). (This function is smooth in the case of a nonstrange attractor and discontinuous 
everywhere in the case of an SNA.) This means that the number of different attractors is equal 
to the number of different iVth iterates, i.e., different ^Ar-values. By setting a tolerance of ±10 -6 
the number of attractors was determined numerically. 

To identify the emergence of SNAs, we compute the attractors and quantify their smooth- 
ness properties using the so-called phase sensitivity exponent introduced by Pikovsky and 
Feudel [1994]. By formally differentiating Eq. (||) with respect to the external phase 1?, we 
get 

dx n+ i „ _ ,dx r . 



(1 + K cos 2irx n ) — — - + 2-ire cos 27r$ n . 
ov 



(6) 



When Eq. (||) is iterated together with Eqs. (§)-( 
and Oxo/dfi = 0, the phase sensitivity 



starting from some initial point (xq,i?o 



Tjv = min max 

(so,i?o) 0<n<N 



dx n 



(7) 
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Figure 21: Minimum number of iterations Fk for which the phase sensitivity saturates, as a function of e, with 
K — 0.9 and Q — 0. Compare with Fig. |(| where k — 33 (dashed line) is used as a threshold. In the yellow marked 
interval the phase sensitivity always saturates for k — 23, indicating two smooth attracting invariant curves. In 
the white marked interval with e smaller than ~ 1.566 saturations occurs for k < 41, indicating a quite wrinkled 
but smooth attracting invariant curve. Using a maximum of F43 iterations, there is no saturation for larger e. 



diverges like for large N in the case of an SNA. On the other hand, in the case of a smooth 
attractor it saturates, i.e. the phase sensitivity exponent \x is zero. The criteria for saturation 
we employ are Tjy < 10 15 and /i < 0.25 (obtained by fitting the slope in a In- In diagram using 
three different iVs) . The black area in Fig. |6| shows the parameter region for which T jy does not 
saturate when eight different initial points are iterated for up to N = F33 = 3 524 578 time steps. 
How much of this black area persists as iV tends to infinity? Figure ^ provides more numerical 
results for larger N and fixed K = 0.9, showing the minimum number of iterations N = (in 
terms of the index of the Fibonacci number) for which the phase sensitivity saturates versus 
the forcing amplitude e. Although a large number of iterations, -F43 = 433494437, was used, 
saturation occurs only for e < 1.566 and e > 1.568. Hence the interval e £ [1.566, 1.568] contains 



the nonsmooth pitchfork bifurcation point, as described in Sec. 4.2, which is approximately 
1.5676. In the next section we discuss the numerical determination of the length of this gap. 



A. 3 Numerical issues regarding the nonsmooth pitchfork bifurcation 

Figure ^ shows the minimum number of iterations N = (in terms of the index of the 
Fibonacci number) for which the phase sensitivity exponent converges, versus e. There is no 
convergence of the phase sensitivity exponent, even for N as large as -F43 = 433494437, in the 
interval e G [1.566, 1.568]. Close to the right side of this interval, a pitchfork bifurcation occurs, 
because there is one smooth attractor (white region in Fig. 21) for smaller e and there are two 
smooth attractors (yellow region) for larger e. Hence, we clearly have a gap when we bound 
the number of iterations by F\~ = F43, but it is unclear whether this gap has nonzero width as 
k — ► 00. 

There are other methods to assess the smoothness properties of the attractor (s). For example, 
in Fig. 22 we applied the method of rational approximations. This method is based on the 
approximation of the irrational frequency uj by rational frequencies tok = Fk—i/Fk with k £ N 
and u) = lmifc^oo Wfc, replacing the quasiperiodically forced map ©-© with a sequence of 
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1.562 1.564 1.566 1.568 1.570 



Figure 22: Minimum order k of the rational approximation for which the sensitive dependence on the initial 
phase #o vanishes (lower curve) and for which the maximum derivative of the attracting sets in #o £ [0, 1) 
saturates (upper curve), as a function of e with K = 0.9 and Q, — 0, cf. Fig. We used 1000 different values of 
0o6[O,l/F fc ). 



periodically (with period forced maps. The F^th iteration of such a map is an orientation- 
preserving diffeomorphism on a circle depending on (£l,£,K) and on the initial phase i?o- The 
union of all attracting invariant sets of this family of diffeomorphisms with $o £ [0, l/F^), forms 
the kth approximation of the attractors of the quasiperiodically forced system. (It is sufficient 
to consider the subinterval [0, 1/F/.), since diffeomorphisms with #o £ [n/F^, (n + 1)/Fk), n = 
1,2,... , Ff; — 1, are topologically conjugate.) For smooth attractors there is a number for which 
the rational approximation of order k and larger does not depend sensitively on i?q [ Pikovsky fc 



Feudel. 



19941. 
1 



Figure [22| shows that this is here the case: the dependence on the initial phase for moderate 
k — one or two stable fixed points for all e, depending on #o — disappears as we cross the 
lower curve towards higher k values — one stable fixed point in the white marked interval, 
corresponding to a single smooth attractor in the quasiperiodically forced map, and two stable 
fixed points in the yellow marked interval, corresponding to a pair of smooth attractors. 

However, as mentioned in pikovsky fc Feudel, 1994 1, a vanishing dependence on i?o is a nec- 
essary but not sufficient condition for smoothness. We also have to stipulate that the maximum 
derivative of the attracting sets with respect to i?o is bounded for all $o £ [0, 1) as k — > oo. 
Note that it is nevertheless sufficient to determine the attracting sets only in the subinterval 
[0, 1/ffc) since the other parts can be obtained by iterating the map ©-(|3|) with the rational 
frequency lou- Furthermore, it is elegant to iterate Eq. (^) simultaneously in order to determine 
the derivative. The maximum derivative obtained by this procedure is an approximation of the 
phase sensitivity. Using the same criteria for boundedness as before, we surprisingly find that 



the maximum derivative always saturates at some order k, shown by the upper curve in Fig. 22 



The qualitative features of Figs. |21|-|22| persist if we choose different Ks in the region with large 
phase sensitivity. 
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A. 4 SNAs in the third region of overlap 

To compute SNAs near the boundary of the tongue in the third region of overlap we use a more 



effective method than the method described in App. A. 2 . This method takes advantage of the 
fact that the corresponding SNAs are unbounded in the x-direction in the lift of the map. We 
define the amplitude of an attractor at time ./V as 

Tn = m in ( max x n — min x r 

(xo,0o) \0<n<N 0<n<AT 



For an unbounded SNA the asymptotic behavior of the amplitude is given by c In A [ Feudel et al., 



1995], in all other cases the amplitude saturates for large N. The numerical algorithm is almost 
the same as for the phase sensitivity (0). However, by fitting the slope using four different N, 
we find as criteria for saturation rV < 6.0 and c < 0.02. Therefore, typically far less iterations 
are necessary, so that a larger maximum number of iterations N = F^q = 102 334 155 and 10 
different initial conditions can be used. 



